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ABSTRACT 

An essential ingredient in kinematic dynamo models of the solar cycle is the internal 
velocity field within the simulation domain - the solar convection zone. In the last 
decade or so, the field of helioseismology has revolutionized our understanding of this 
velocity field. In particular, the internal differential rotation of the Sun is now fairly 
well constrained by helioseismic observations almost throughout the solar convection 
zone. Helioseismology also gives us some information about the depth-dependence of 
the meridional circulation in the near-surface layers of the Sun. The typical velocity 
inputs used in solar dynamo models, however, continue to be an analytic fit to the ob- 
served differential rotation profile and a theoretically constructed meridional circulation 
profile that is made to match the flow speed only at the solar surface. Here we take the 
first steps towards the use of more accurate velocity fields in solar dynamo models by 
presenting methodologies for constructing differential rotation and meridional circula- 
tion profiles that more closely conform to the best observational constraints currently 
available. We also present kinematic dynamo simulations driven by direct helioseismic 
measurements for the rotation and four plausible profiles for the internal meridional 
circulation - all of which are made to match the helioseismically inferred near-surface 
depth-dependence, but whose magnitudes are made to vary. We discuss how the results 
from these dynamo simulations compare with those that are driven by purely analytic 
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fits to the velocity field. Our results and analysis indicate that the latitudinal shear in 
the rotation in the bulk of the solar convection zone plays a more important role, than 
either the tachocline or surface radial shear, in the induction of the toroidal field. We 
also find that it is the speed of the equatorward counterflow in the meridional circula- 
tion right at the base of the solar convection zone, and not how far into the radiative 
interior it penetrates, that primarily determines the dynamo cycle period. Improved 
helioseismic constraints are expected to be available from future space missions such 
as the Solar Dynamics Observatory and through analysis of more long-term continuous 
data sets from ground-based instruments such as the Global Oscillation Network Group. 
Our analysis lays the basis for the assimilation of these helioseismic data within dynamo 
models to make future solar cycle simulations more realistic. 

Subject headings: Sun: magnetic fields. Sun: interior. Sun: rotation. Sun: helioseismol- 
ogy. Sun: activity 

1. Introduction 

The dynamic nature of solar activity can often be traced back to the presence and evolution 
of magnetic fields in the Sun. The more intense magnetic fields on the order of 1000 Gauss (G) 
are observed to be concentrated within regions known as sunspots, which often appear in pairs of 
opposite magnetic polarities (Hale 1908). Sunspots have been observed regularly now for about 
four centuries starting with the telescopic observations of Galileo Galilei in the early 1600s. These 
observations point out that the number of sunspots on the solar surface varies in a cyclic fashion 
with an average periodicity of 11 years (Schwabe 1844), although there are variations both in the 
amplitude and period of this cycle. At the beginning of a cycle sunspots appear at about mid- 
latitudes in both hemispheres (with opposite polarity orientation across the hemispheres) and then 
progressively appear at lower and lower latitudes as the cycle progresses (Carrington 1858) until no 
sunspots are seen (i.e., solar minimum). In the next cycle, the same pattern repeats, but the new 
cycle spots have their bipolar magnetic orientation reversed relative to the previous cycle (in both 
hemispheres). So considering sign as well as amplitude, the solar cycle has a period of 22 years. 

There is a weaker, more diffuse component of the magnetic field outside of sunspots which 
is seen to have a somewhat different evolution. This field - whose radial component has been 
observable at the solar surface since the advent of the magnetograph - was believed to be on the 
order of 10 Gauss. It is found that this field is concentrated in unipolar patches, which originate 
at sunspot latitudes at the time of sunspot maxima, and then moves poleward with the progress 
of the sunspot cycle (Babcock 1959; Bumba and Howard 1965; Howard and LaBonte 1981). The 
sign of this radial field of any given cycle is opposite to the old cycle polar field, which it cancels 
and reverses upon reaching the poles. The amplitude of this radial field achieves a maximum 
(at the poles) at the time of sunspot minimum (i.e. with a 90° phase difference relative to the 
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sunspots). However, the periodicity of the cycle of this field matches the sunspot cycle period, 
underscoring that they are related. Recent observations by Hinode indicate that this radial field 
gets concentrated within unipolar flux tubes with field strength on the order of 10^ Gauss (Tsuneta 
et al. 2008). 

Explanations of this observations of the solar magnetic cycle rely on the field of magneto- 
hydrodynamic dynamo theory, which seeks to address the generation and evolution of magnetic 
fields as a complex non-linear process involving interactions between the magnetic field and plasma 
flows within the solar interior. In particular, it is now believed that the solar cycle involves the 
generation and recycling (feeding on the energy available in plasma motions) of two components 
of the magnetic field - the toroidal component and the poloidal component. In an axisymmetric 
spherical coordinate system, the magnetic and velocity fields can be expressed as 



where the first term on the right hand side (R.H.S.) of Equations 1 and 2 is the toroidal component 
(in the case of the velocity field this corresponds to the differential rotation) and the second term 
is the poloidal component of the field (in the case of the velocity field this corresponds to the 
meridional circulation). The toroidal component of the magnetic field is thought to be produced 
by stretching of an initially poloidal field by the differential rotation of the Sun (the dynamo 
ri-effect); subsequently, strong toroidal flux loops rise up due to magnetic buoyancy emerging 
through the solar surface as sunspots (Parker 1955a). To complete the dynamo cycle, the poloidal 
field (whose radial component is manifested as the observed vertical field on the solar surface) 
has to be regenerated from this toroidal field in a process that is traditionally called the dynamo 
a-effect. The first explanation of this a-effect was due to Parker (1955b) who suggested that helical 
turbulent convection in the solar convection zone (SCZ) would twist rising toroidal flux tubes into 
the poloidal plane, recreating the poloidal component of the magnetic field. Much has changed 
since this pioneering description of the first solar dynamo model by Parker, although the basic 
notion of the recycling of the toroidal and poloidal components remain the same. 

First of all, simulations of the buoyant rise of toroidal flux tubes point out that to match 
the observed properties of sunspots at the solar surface, the strength of these flux tubes at the 
base of the SCZ has to be much more than the equipartition field strength of 10^ G (D'Silva 
& Choudhuri 1993; Fan, Fisher & DeLuca 1993). The classical dynamo a-effect due to helical 
turbulence is expected to be quenched for super-equipartition field strengths and therefore other 
physical processes have to be invoked as a regeneration mechanism for the poloidal field. One of 
the alternatives is an idea originally due to Babcock (1961) and Leighton (1969). The Babcock and 
Leighton (hereby BL) model proposes that the decay and dispersal of tilted bipolar sunspot pairs at 
the near-surface layers, mediated by diffusion, differential rotation, and meridional circulation, can 
regenerate the poloidal field. This mechanism is actually observed and is simulated as a surface flux 
transport process that can reproduce the solar polar field reversals (Wang, Nash & Sheeley 1989). 



B = B^e^ + V X (Ae^) 
V = r sin(^)ile(^ + Vp 



(1) 

(2) 
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Therefore a synthesis of Parker's original description along with the BL mechanism for poloidal 
field generation is now widely accepted as a leading contender for explaining the solar dynamo 
mechanism (Choudhuri, Schiissler &; Dikpati 1995; Durney 1997; Dikpati &; Charbonneau 1999; 
Nandy & Choudhuri 2001; Nandy 2003), although there are other alternative suggestions as well. 
A description of all of those is beyond the scope of this paper and interested readers are referred 
to the review by Charbonneau (2005). 

Second, helioseismology has now mapped the solar internal rotation profile (Schou et al. 1998; 
Charbonneau et al. 1999), which is observed to be vary mainly in the latitudinal direction in the 
main body of the SCZ. Helioseismology has also discovered the tachocline - a region of strong radial 
and latitudinal shear beneath the base of the SCZ which is expected to play an important role in 
the generation and storage of strong toroidal flux tubes. 

Third, more is now known about the meridional circulation, which is observed to be polewards 
at the surface (Hathaway 1996). To conserve mass, this circulation should turn equatorwards in 
the solar interior. This circulation is deemed to be important for the dynamics of the solar cycle 
(see Hathaway et al. 2003 and the review by Nandy 2004) but the profile of this in the solar 
interior remains poorly constrained. Nandy and Choudhuri (2002) proposed a deep equatorward 
counterflow in the circulation (penetrating into the radiative interior beneath the SCZ) to better 
reproduce in dynamo simulations the latitudinal distribution of sunspots (because equatorward 
advective transport and storage of the deep seated toroidal field is more efficient at these depths 
where turbulence is greatly reduced). However, Gilman and Miesch (2004), based on a laminar 
analysis, argued that the penetration of the circulation would be limited to a shallow Ekman layer 
close to the base of SCZ. A recent and more detailed analysis of the problem by Garaud and 
Brummell (2008) suggests that the circulation can in fact penetrate deeper down into the radiative 
interior. At this point there is no consensus on the profile and nature of the meridional circulation 
in the solar interior. Helioseismic data does provide some information about the depth-dependence 
of this circulation at near-surface layers (Braun & Fan 1998, Giles 2000, Chou & Ladenkov 2005 , 
Gonzalez-Hernandez et al. 2006), which, in conjunction with reasonable theoretical arguments, can 
be used to construct some plausible interior profiles of this flow. 

Numerous kinematic dynamo models have been constructed in recent years (see Charbonneau 
2005 for a review) incorporating these large scale flows (differential rotation and meridional cir- 
culation) as drivers of the magnetic evolution. More recently such dynamo models (based on the 
BL idea of poloidal field generation) have also been utilized to make predictions for the upcoming 
cycle (Dikpati, de Toma & Gilman 2006; Choudhuri, Chatterjee & Jiang 2007). At present, all 
these kinematic dynamo models incorporate the information on large-scale flows as analytic fits to 
the differential rotation profile and a theoretically constructed meridional circulation profile that 
is subject to mass conservation but matches the flow speed only at the solar surface (i.e., without 
incorporating the depth-dependent information that is available). However, these large-scale flows 
are crucial to the generation and transport of magnetic fields; the differential rotation is the primary 
source of the toroidal field that creates solar active regions, and the meridional flow is thought to 
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play a crucial role in coupling the two source regions for the poloidal and toroidal field through 
advective flux transport. Given this, it is obvious that the next step in constructing more sophis- 
ticated dynamo models of the solar cycle is to move towards a more rigorous use of helioseismic 
data to constrain these models in a way such that they conform more closely to the best available 
observational constraints; that is the goal of this study 

In Section 2, we describe the basic features of the kinematic dynamo model based on the BL 
idea that we use for our study; in this model, we use fairly standard parametrization (commonly 
used in the community) of various processes such as the diffusivity, dynamo a-effect and magnetic 
buoyancy. In Sections 3.1 and 3.2 we present the methodologies for using the helioseismically 
observed solar differential rotation and constraining the meridional circulation profiles within this 
dynamo model and describe how they improve upon the commonly used analytic profiles. In Section 
4 we present results from dynamo simulations using these improved helioseismic constraints and 
conclude in Section 5 with a summary of our main results and their contextual relevance. 



Our Model 



We substitute Equations 1 and 2 into the magnetic induction equation 

— = V X (v X B - ryV X B) 

dt ^ ' ^ 



(3) 



and add the phenomenological BL poloidal field source a to obtain the axisymmetric dynamo 
equations: 



dA 



dB 

~dt 



+ s 



Vp- V 



+ {V'Vp)B = r/ V 



B + s{[\/ X {Ae^)]'\/n) 



(4) 



(5) 



1 djsB) dr] 
s dr 9r ' 



where s = rsin(^). The terms on the LHS of both equations with the poloidal velocity (vp) 
correspond to the advection and deformation of the magnetic field by the meridional flow. The first 
term on the RHS of both equations corresponds to the diffusion of the magnetic field. The second 
term on the RHS of both equations is the source of that type of magnetic field (BL mechanism for 
A and rotational shear for B). Finally, the third term on the RHS of Equation 5 corresponds to 
the advection of toroidal field due to a turbulent diffusivity gradient. 

As mentioned before, the generation of poloidal field near the surface due to the decay of active 
regions is modeled through the inclusion of a source term, which acts as a source for the vector 
potential A. It is localized both in radius and latitude matching observations of active region 
emergence patterns and depends on the average toroidal field at the bottom of the convection 
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zone (Dikpati & Charbonneau 1999). The radial and latitudinal dependence of the source is the 
following: 



0) = co,(«) (l + erf (O-i^O'-^A ) (, - erf " <»° + 



16 \ \ 7 //V V 7 

dal J J \ V "a/i 



(6) 



where q;o sets the strength of the source term and we set it to the value in which our system 
starts having oscillating solutions. The parameters f3 = 40^ and 7 = 10^ characterize the latitudes 
in which sunspots appear. On the other hand, r^/ = O.94i?0, dai = O.O4i?0, Tah — and 
dah — O.Oli?0 characterize the radial extent of the region in which the poloidal field is deposited 
(see Figure [1]- a and b). Besides radial and latitudinal dependence, we also introduce lower and 
upper operating thresholds on the poloidal source that is dependent on toroidal field amplitude 
- which is a more realistic representation of the physics involving magnetic buoyancy (as argued 
in Nandy & Choudhuri 2001; Nandy 2002; see also Charbonneau, St- Jean & Zacharias 2005). 
The presence of a lower threshold is in response to the fact that the plasma density inside weak 
flux tubes is not low enough (compared to the density of the surrounding plasma) to make them 
unstable to buoyancy (Caligari, Moreno-Insertis & Schussler 1995) and those that manage to rise 
have very long rising times (Fan, Fisher and De Luca 1993). On the other hand, if flux tubes are 
too strong they are not tilted enough when they reach the surface to contribute to poloidal field 
generation (D'Silva and Choudhuri 1993; Fan, Fisher and Deluca 1993). The dependence of the 
poloidal source on magnetic field is 

FiBa,) = ^ 2 ( 1 ^ > (7) 

where Kae = 1/ niax(F(Sa-y)) is a normalization constant and = 1.5 x lO^G and Bi — 4:X lO^G 
are the operating thresholds (see Figure [l]-c) . 

Another ingredient of this model is a radially dependent magnetic diffusivity; in this work we 
use a double-step profile (see Figure [2]) given by 

Vir) = Vbcd + ^1 + erf J J + ^1 + erf , J J (8) 

where rii,cd — W^crri^ /s corresponds to the diffusivity at the bottom of the computational domain, 77C2 - 
lO-*^-*^ corresponds to the diffusivity in the convection zone, rjsg = 10^^cm?/s corresponds to 
the supergranular diffusivity and Tcz = O.73i?0, dcz = O.O3i?0, Tgg = O.95i?0 and dgg = O.O5i?0 
characterize the transitions from one value of diffusivity to the other. Although it is common these 
days to use these sort of profile, a point we note in passing is that the exact depth dependence 
of turbulent diffusivity within the SCZ is poorly, if at all, constrained. Besides the value of the 
supergranular diffusivity r]sg, which can be estimated by observations (Hathaway & Choudhary 
2008), the others parameters are not necessarily well constrained. 
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Once all ingredients are defined(see Section|3]for the flow field), we solve the dynamo equations 
using a recently developed and novel numerical technique called exponential propagation (see Ap- 
pendix). Our computational domain is defined in a 250 x 250 grid comprising only one hemisphere. 
Since we run our simulations only in one hemisphere our latitudinal boundary conditions at the 
equator {9 — n /2) are dA/dO — and B — Q. Furthermore, since the equations we are solving are 
axisymmetric, both the potential vector and the toroidal field need to be zero (A = and S = 0) at 
the pole (0 = 0), to avoid singularities in spherical coordinates. For the lower boundary condition 
(r = O.55i?0), we assume a perfectly conducting core, such that both the poloidal field and the 
toroidal field vanish there (i.e., A, S = at the lower boundary). For the upper boundary condition 
we assume that the magnetic field has only a radial component {B — Q and d{rA)/dr — 0); this con- 
dition has been found necessary for stress balance between subsurface and coronal magnetic fields 
(for more details refer to van Ballegooijen and Mackay 2007). As initial conditions we set A = 
throughput our computational domain and B oc sin(20) * sin(7r * ((r — O.55i?0)/(i?0 — O.55i?0))). 
After a few cycles, all transients related to the initial conditions typically disappear and the dynamo 
settles into regular oscillatory solutions whose properties are determined by the parameters in the 
dynamo equations. 



3. Constraining the Flow Fields 

The last two ingredients of the dynamo model are the velocity fields (differential rotation and 
meridional circulation), which are the focus of this work. The differential rotation is probably the 
best constrained of all dynamo ingredients but the actual helioseismology data has never before 
been used directly in dynamo model, only an analytical fit to it. We discuss below how the 
actual rotation data can be used directly within dynamo models through the use of a weighting 
function to filter out the observational data in the region where it cannot be trusted. On the other 
hand, the meridional circulation is one of the most loosely constrained ingredients of the dynamo. 
Traditionally only the peak surface flow speed is used to constrain the analytical functions that are 
used to parameterize it, in conjunction with mass conservation. In this work we take advantage 
of the properties of such functions and make a fit to the helioseismic data on the meridional flow 
that constrains the location and extent of the polar downflow and equatorial upflow, as well as the 
radial dependence of the meridional flow near the surface - thereby taking steps towards better 
constrained flow profiles. 



3.1. Differential Rotation 

As opposed to meridional circulation, there are helioseismic measurements of the differential 
rotation for most of the convective envelope which can be used directly in our simulations. Here 
we use data from the Global Oscillation Network Group (GONG) (courtesy Dr. Rachel Howe) 
obtained using the RLS inversion mapped onto a 51 x 51 grid (see Figure Isla). However, these 
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observations cannot be trusted fully in the region within O.3i?0 of the rotation axis (specifically at 
high latitudes), because the inversion kernels have very little amplitude there. Below we outline a 
method to deal with this suspect data by creating a composite rotation profile that replaces these 
data at high latitudes with plausible synthetic data, that smoothly matches to the observations in 
the region of trust. 



3.1.1. Adaptation of the Data to the Model 

In the first step we use a splines interpolation in order to map the data to the resolution of 
our simulation (a grid of 250 x 250 see Figure [sj-a). The next step is to make a composite with the 
data and the analytical form of Charbonneau et al. (1999; see Figure [sj-b). The analytical form is 
defined as: 



f^A(r, e) = 2^ [n, + ^ (l - erf (^)) {^e - f^c + - ^e)^s{0)) 



(9) 



^s{0) = acos^{e) + (1 - a) cos4(0), 

where — 432 nHz is the rotation frequency of the core, fig = 470 nHz is the rotation frequency 
of the equator, — 330 nHz is the rotation frequency of the pole, a — 0.483 is the strength of the 
cos^(0) with respect to the cos^(0) term, rtc = 0.716 the location of the tachocline and wtc = 0.03. 
We use the parameters defining the tachocline 's location and thickness as reported by Charbonneau 
et al. (1999) for a latitude of 60^. This is because they match the data better at high latitudes 
(which is the place where the data merges with the analytical profile) than those reported for the 
equator. This composite replaces the suspect data at high latitudes within O.3i?0 of the rotation 
axis with that of the analytic profile. However, it is important to note that at low latitudes, within 
the convection zone, the actual helioseismic data is utilized. 

In order to make the composite we create a weighting function m(r, 9) with values between 
and 1 for each grid-point defining how much information will come from the RLS data and how 
much from the analytical form (see Figure [sj-d) . We define the weighting function in the following 
way: 



m(r,^) = 1- M 1 + erf ' ^^^^T^ , (10) 



1 / /r2cos(2^) -r2 

where Vm — O.5i?0 is a parameter that controls the center of the transition and dm — O.6i?0 controls 
the thickness. The resultant differential rotation profile, which can be seen in Figure [3]-c, is then 
calculated using the following expression: 



Q(r, e) = m{r, 9)nRLs{r, 9) + [I - m{r, 9)]nA{r, 6) 



(11) 
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3.1.2. Differences Between the Analytical Profile and the Composite Data 

It is instructive to compare the analytical and composite profiles with the actual helioseis- 
mology data. In Figure [4]-a we present the residual error of subtracting the analytical profile of 
Charbonneau et al. (1999), with rtc = 0.7 and wtc = 0.025, from the RLS data. In Figure [4]-b we 
present the residual error of subtracting our composite from the RLS data. As is expected, there 
is no difference between the composite and raw data at low latitude, but the residual increases as 
we approach the rotation axis - where the RLS data cannot be trusted. For the analytical profile, 
the residuals errors are more significant, even at low latitudes. This demonstrates the ability of our 
methodology to usefully integrate the helioseismic data for differential rotation. 



3.2. Meridional Circulation 

The meridional flow profile remains rather poorly constrained in the solar interior even though 
the available helioseismic data can be used to constrain the analytic flow profiles that are in use 
currently. Here we present ways to betters constrain this profile with helioseismic data. In order 
to do that, we use data from GONG (Courtesy Dr. Irene Gonzalez-Hernandez) obtained using the 
ring-diagrams technique. This data, which we can see in Figure |6]-a, corresponds to a time average 
of the meridional flow between 2001 and 2006; and it comprises 19 values of r from Rq down to 
a depth of O.97i?0, and 15 different latitudes between —52.5^ and 52.5^. It is important to note 
that our work relies heavily in the assumption that the meridional flow is adequately described 
by a stream function with separable variables. This is consistent with the assumption present 
implicitly in all work on axisymmetric solar dynamo models up to this date. Below we use this 
property of our stream function, along with weighted latitudinal and radial averages of the data, 
to completely constrain its latitudinal dependence, as well as the topmost ten percent of its radial 
dependence. As this data currently does not constrain the depth of penetration of the flow in the 
deep solar interior, we explore two different plausible penetration depths of the circulation. For 
reasons described later, we choose to perform simulations with two different peak meridional flow 
speeds, therefore exploring four plausible meridional flow profiles altogether. 



3.2.1. Constraining the Latitudinal Dependence of the Meridional Flow 

Meridional circulation has been typically implemented in these type of dynamo models by 
using a stream function in combination with mass conservation, i.e.: 

Vp(r,e) = -^Vx(v|>(r,e)e^). (12) 

The two stream functions that are commonly used were proposed by van Ballegooijen and Choud- 
huri (1988) and by Dikpati and Choudhuri (1995). They have in common the separability of 
variables and thus can be written in the following way: 
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^{r,e)^voF{r)G{e), (13) 
where 'Uq is a constant which controls the amphtude of the meridional flow. 

Using such a stream function the components of the meridional flow become: 

"»(''.<') = -"«;:^|:(''-^M)G«'). (15) 

which can themselves be separated into the multiplication of exclusively radially and latitudinally 
dependent functions. This property allows us to constrain the entire latitudinal dependence of this 
family of functions by using the available helioseismology data for vq at the surface. This can be 
done because the latitudinal dependence of vq is exactly the same as that of the stream function 
and only the amplitude of this functional form changes with depth, see for example Figure [5] for 
the latitudinal velocity at different depths used by van Ballegooijen and Choudhuri (1988). In this 
work we assume a latitudinal dependence like the one they used, i.e., 

G{e) =sin(^+i)(0)cos(0). (16) 

In order to estimate the parameter q we first take a density weighted average of the helioseismic 
data using the values for solar density from the Solar Model S (Christensen-Dalsgaard et al. 1996), 
such that 

MO,) - ^^,(,^) • (17) 

In Figure [6]-b we plot the meridional flow at different depths weighted by density, which we add for 
each latitude in order to find the average. We then use this average to make a least squares fit to 
the analytical expression (eq 16), which we can see in Figure [6[c. We find that a value of g = 1 fits 



the data best. This therefore constrains the latitudinal (6) dependence of the flow profile. 



3.2.2. Constraining the Radial Dependence of the Meridional Flow 

As opposed to the latitudinal dependence, the radial dependence of the meridional flow is less 
constrained since there is no data below O.97i?0. However, at least some of the parameters can be 
constrained: We start with the solar density, for which we perform a least squares fit to the Solar 
Model S using the following expression: 



r 



p(r)~(^-7) (18) 



we find that values of 7 = 0.9665 and m = 1.911 fit the model best (see Figure [7[c). 
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In the second step we constrain the radial dependence of the stream function. We begin with 
the function 



where Rq corresponds to the solar radius, Rp to the maximum penetration depth of the meridional 
flow and a and Ri control the location in radius of the poleward peak and the value of the meridional 
flow at the surface. In order to constrain them we use the helioseismic data again, but this time 
we use the radial dependence (see Figure [7]-a). We flrst remove the latitudinal dependence, which 
we do by dividing the data of each latitude by G{9) using the value of g = 1 found in the previous 
section. In Figure [7]-b we plot the flow data after removing the latitudinal dependence, note that 
there is no longer any sign difference between the two hemispheres. The next step is to generate the 
latitudinal average which we can see as black dots in Figure [7|-c. It is evident from looking at the 
radial dependence of the meridional flow that the velocity increases with depth for most latitudes, 
and that the point of maximum velocity is not within the depth up to which the data extends. 
Since the exact radial dependence of the data is too complex for our functional form to grasp, 
the features we concentrate on reproducing are the presence of a maximum inside the convection 
zone, as well as the amplitude of the flow at the near surface-layers. The logic here is to use the 
fewest possible parameters and a simple, physically transparent proflle that does a reasonable job 
of matching the data. In view of the lack of better constraints, we assume here that the peak of the 
return flow is at O.97i?0 (which is the depth at which the current helioseismic data has it's peak). 

Following the procedures and steps above we construct proflles to flt the depth-dependence 
pointed out in the available helioseismic data; however, this does not constrain how far the merid- 
ional flow can penetrate and therefore we try two different penetration depths, one shallow O.Tli?© 
(i.e., barely beneath the base of the SCZ) and one deep Rp = O.64i?0 (into the radiative interior) - 
both of which match the latitudinal and radial constraints as deduced from the near-surface helio- 
seismic data. Note that observed light element abundance ratios limit the depth of penetration of 
the circulation to about Rp = O.62i?0 (Charbonneau 2007). Now we know that the meridional flow 
speed is highly variable, with fluctuations that can be quite signiflcant and the measured flow speed 
can change depending on the phase of the solar cycle (Hathaway 1996; Gizon & Rempel 2008). The 
magnetic flelds are also expected to feed back on the flow (Rempel 2006). Taken together these 
considerations point out that the effective meridional flow speed to be used in dynamo simulations 
could be less than that implied from the Gonzalez-Hernandez et al. data, a possibility that is borne 
out by the work of Braun & Fan (1998) and Gizon & Rempel (2008), who flnd much lower peak 
flow speeds in the range 12-15 m/s. Keeping this in mind, we use the same latitudinal and radial 
constraint as deduced earlier, but consider in our simulations two additional proflles with peak flow 
speeds of 12 m/s with deep and shallow penetrations. Therefore, in total we explore four plausible 
meridional flow proflles in our dynamo simulations (see Figure [7]-c and Table [l] for an overview) , 
the results of which are presented in the next section. 




(19) 
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Rp 


a 




Set 1 


12 m/s 


OMRq 


1.795 


l.O27i?0 


Set 2 


22 m/s 


Set 3 


12 m/s 


O.IIRq 


2.03 


l.O3i?0 


Set 4 


22 m/s 



Table 1: Sets of parameters characterizing the different meridional flow profiles used in our dynamo 
simulations. Vq corresponds to the meridional flow peak speed, Rp the maximum penetration of 
the flow, and a and Ri are parameters that control the location of the poleward flow as well as the 
surface speed. 



4. Dynamo Simulation: Results and Discussions 
4.1. Analytic vs. Helioseismic Composite Differential Rotation 

We first compare dynamo solutions found using the helioseismology composite of the differential 



rotation (Section |3.1.1| ) and the analytical profile of Charbonneau et al. (1999). In doing so, 
we perform simulations with a model as described in Section 2 (see also numerical methods in 
the Appendix) and generate field evolution maps for the toroidal and poloidal fields. From the 
simulated butterfly diagrams for the evolution of the toroidal field at the base of the SCZ and 
the surface-radial field evolution (Figure [s]), we find that large-scale features of the simulated solar 
cycle are generally similar across the two different rotation profiles (even with different meridional 
flow penetration depths and speeds), specially for the shallowest penetration (sets 3 and 4 with 
Rp = O.71i?0). 



In order to understand this similarity, it is useful to look at Figures [T0| and [TT] corresponding to 
simulations using the composite differential rotation, and the meridional flow sets 1 (Rp = O.64i?0, 
Vo = 12m/s) and 4 {Rp = O.Tli?©, Vq = 22m/s). The first two columns, from left to right of 
both figures show the evolution of the shear sources (B^ • Vr-f^) and (B^ • V^f^) - which contribute 
to toroidal field generation by stretching of the poloidal field. It is evident that the location and 
strength of these sources is different - the radial shear source is mainly present near the surface 
whereas the latitudinal shear is spread throughout the convection zone. It can also be seen that the 
radial shear source is roughly five times stronger than the latitudinal. However, if attention is paid 
to the evolution of the toroidal field (third column from the right in both figures), it is clear that 
this radial shear term has no significant impact on the structure and magnitude of the toroidal field 
(a similar result was found by Dikpati et al. 2002). The reason is that the upper boundary condition 
(S = at r = ^0), in combination with the high turbulent diffusivity (and thus short diffusive 
time scale) there, imposes itself very quickly on the toroidal field generated by the radial shear - 
washing it out. This greatly reduces the relative role of the surface shear as a source of toroidal 
magnetic field, effectively making the surface dynamics very similar across simulations using the 
analytical rotation profile (without any surface shear) or the composite helioseismic profile. 
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Once the surface layers are ruled out as important sources of toroidal field generation we are 
faced with the fact that the strongest source of toroidal field is the latitudinal shear inside the 



convection zone (and not the tachocline radial shear), as is evident in Figures 10 and 11 This 
goes against the commonly held perception that the tachocline is where most of the toroidal field is 
produced. However, the importance of the latitudinal shear term in the SCZ is clearly demonstrated 
here where we have plotted the shear source terms, which is not normally done by dynamo modelers. 

The establishment of the SCZ as an importance source region of toroidal field is of relevance 
when regarding the similarity of the solutions obtained using the composite data and the analytical 
profile. This is because the region where the shear of the analytical profile and the helioseismic 
composite data differ most (both for radial and latitudinal shear) is in the tachocline. This is 
evident in Figure [9] where we plot the residual of subtracting the radial and latitudinal shear of 
the analytical profile from the shear of the composite data. The similarity between solutions is 
specially important for shallow meridional flow profiles with low penetration - which does not 
transport any poloidal field into the deeper tachocline, thereby further diminishing the role of the 
tachocline shear. 



4.2. Shallow versus Deep Penetration of the Meridional Flow 

In the second part of our work we compare dynamo solutions obtained for each of the four 
different meridional flow profiles with two different penetration depth and with two different peak 
flow speeds. First, from Figure [8] it is evident that the shape of the solutions changes with varying 
penetration depth; this is caused by the increasing role of the tachocline shear in generation and 
storage of the toroidal field as the penetration depth of the meridional flow increases. This is 
apparent when one compares Figure [To] (deep penetration) to Figure pT] (shallow penetration). If 
we look at the inductive shear sources it is evident that for the shallowest penetration no field is 
being generated inside the tachocline. In the poloidal field plots (right column) we see that no 
poloidal field is advected into the tachocline region for the shallowest penetration, but some is 
advected for the deepest. 

Second we compare the periods of our solutions in Table [2j We find that most solutions have 
a sunspot cycle (i.e., half-dynamo cycle) period that is comparable to that of the Sun, with the 
exception of the fast flow with deep penetration and the slow flow with shallow penetration, which 
have respectively a comparatively smaller and larger period. As the meridional flow is buried 
deeper, one expects the length of the advective circuit to increase, thereby resulting in larger 
dynamo periods. However, it is evident that as we increase the penetration, the period decreases 
- even if the length of the flow loop that supposedly transports magnetic flux increases; this is 
counterintuitive but has a simple explanation. Our simulations and exhaustive analysis points out 
that it is not how deep the flow penetrates that governs the cycle period, but it is the magnitude of 
the meridional counterflow right at the base of the SCZ (Rp — OJISRq) that is most relevant. This 
is because most of the poloidal field creation at near-surface layers is coupled to buoyant eruptions 
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of toroidal field from this layer of equatorward migrating toroidal field belt at the base of the SCZ 



(see third columns in Figures [10| and 11) and it is the flow speed at this region that governs the 
dynamo period. In Figure [7| it is clear that the speed of the counterflow in this convection zone- 
radiative interior interface increases as the flow becomes more penetrating (a consequence of the 
constraints set by mass conservation and the fits to the near-surface helioseismic data), thereby 
reducing the dynamo period. 

Overall, an evaluation of the butterfly diagram (Figure [s]), points out that the toroidal field 
belt extends to lower latitudes (where sunspots are observed) for deeper penetrating meridional 
flow, although there is a polar branch as well. For the shallow flow, we find that the toroidal field 
belt is concentrated around mid-latitudes with almost symmetrical polar and equatorial branches 
- a signature of the convection zone latitudinal shear producing most of the toroidal field (as in 
the interface dynamo models, see, e.g., Parker 1993 and Charbonneau & MacGregor 1997). 





Rp 


Vo 


r - HS data (yrs) 


r - Analytical (yrs) 


Set 1 


OMRq 


12 m/s 


9.67 


10.00 


Set 2 


22 m/s 


5.63 


5.67 


Set 3 


O.71i?0 


12 m/s 


14.67 


14.02 


Set 4 


22 m/s 


11.85 


12.85 



Table 2: Simulated sunspot cycle period for the different sets of meridional flow parameters. Rp 
corresponds to the maximum penetration depth of the meridional flow, Vq to the peak speed in the 
poleward flow and r is the period of the solutions in units of years. For the rest of the parameters 
in each set please refer to Table [l] 



4.3. Dependence of the solutions on changes in the turbulent diffusivity profile 

Although a detailed exploration of the turbulent diffusivity parameter space is outside the 
scope of this work, we study two special cases in which we vary a single parameter while leaving 
the rest fixed. 

For the first case we lower the diffusivity in the convection zone, tJcz, from lO^^cm^/s to 
lO^^cm^/s. As can be seen in Figure 12, this introduces two important changes in the dynamo 



solutions: The first one is an overall increase in magnetic field magnitude due to the reduction in 
diffusive decay while keeping the strength of the field sources constant. The second is a drastic 
increase in the dynamo period (which can be seen tabulated in Table [s]) for the solutions that use a 
meridional flow with a penetration of Rp = O.71i?0. The reason behind such a change resides in the 
nature of the transport processes at the bottom of the convection zone, which are a combination of 
both advection and diffusion. In the case of flow profiles with deep penetration the velocity at the 
bottom is high enough for downward advection to transport flux into the tachocline. On the other 
hand, in the case of low penetration, the last bit of downward transport into the tachocline is done 
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by diffusive transport and thus dominated by diffusive timescales. Because of this, by decreasing 
turbulent diffusivity by an order of magnitude, we drastically increase the period of the solutions. 





Rp 


Vo 


r - HS data (yrs) 


r - Analytical (yrs) 


Set 1 


OMRq 


12 m/s 


12.46 


12.86 


Set 2 


22 m/s 


6.47 


6.55 


Set 3 


O.71i?0 


12 m/s 


87.78 


90.92 


Set 4 


22 m/s 


70.48 


74.41 



Table 3: Simulated sunspot cycle period for the different sets of meridional flow parameters when 
using a low diffusivity in the convection zone {ijcz = lO^^cm^/s). Rp corresponds to the maximum 
penetration depth of the meridional flow, Vq to the peak speed in the poleward flow and r is the 
period of the solutions in units of years. 

In the second parameter space experiment, we lower the super- granular diffusivity ijsg from 
10 cm /s to IQi^cm^/s. This was done to study the impact of the surface radial shear under low 



diffusivity conditions. However, as can be seen in Figure [131 there is very little difference between 
the two solutions. This means that even after reducing the super-granular diffusivity by two orders 
of magnitude, the radial shear has very little impact on the solutions and the upper boundary 
conditions still play an important role in limiting the relative contribution from the near-surface 
shear layer. 



5. Conclusions 

In summary, we have presented here methods which can be used to better integrate helioseismic 
data into kinematic dynamo models. In particular, we have demonstrated that using a composite 
between helioseismic data and an analytical profile for the differential rotation, we can directly use 
the helioseismic rotation data in the region of trust and substitute the suspect data by smoothly 
matching it to the analytical profile where the data is noisy. This paves the way for including the 
helioseismically inferred rotation profile directly in dynamo simulations. We have also shown how 
mathematical properties of the commonly used analytic stream functions describing the meridional 
flow can be fit to the available near-surface helioseismic data to entirely constrain the latitudinal 
dependence of the meridional flow, as well as weakly constrain the radial (depth) dependence. 

In our simulations, comparing the helioseismic data for the differential rotation with the an- 
alytical profile of Charbonneau et al. (1999), with four plausible meridional flow profiles, we find 
that there is little difference between the solutions using the helioseismic composite and the ana- 
lytical differential rotation profile - specially for shallow penetrations of the meridional flow and 
even at reduced super-granular diffusivity. This is because the impact of the surface radial shear, 
which is present in the helioseismic composite but not the analytic profile, is greatly reduced by 
the proximity of the upper boundary conditions. Also, for the shallow circulation, the toroidal field 



-16- 



generation occurs in a region located above the tachocline with mainly latitudinal shear, where the 
difference between the composite data and the analytical profile is not significant. 

The main result from this comparative analysis is that the latitudinal shear in the rotation is 
the most dominant source of toroidal field generation in these type of models that are characterized 
by high diffusivity at near surface layers, but lower diffusivity within the bulk of the SCZ - specially 
near the base where most of the toroidal field is being created. Since this latitudinal shear exists 
throughout the convection zone, an interesting question is whether toroidal fields can be stored there 
long enough to be amplified to high values by the shear in the rotation, without being removed by 
magnetic buoyancy. If this were to be the case, i.e., the latitudinal shear is indeed confirmed to be 
the dominant source of toroidal field induction, we anticipate then that downward flux pumping 
(Tobias et al. 2001; see also Guerrero & de Gouveia Dal Pino 2008) - which tends to act against 
buoyant removal of flux, may have an important role to play in this context. This could also call 
into question the widely held view that the solar tachocline is where most of the toroidal field 
is created and stored (see Brandenburg 2005 for arguments favoring a more distributed dynamo 
action throughout the SCZ). 

Our attempts to integrate helioseismic meridional flow data into dynamo models and related 
simulations have uncovered points that are both encouraging and discouraging. 

On the discouraging side, we find that the currently available observational data are inade- 
quate to constrain the nature and exact profile of the deep meridional flow, especially the return 
flow. Neither do the simulation results and their comparison with observed features of the solar 
cycle clearly support or rule out any possibility. A recent analysis on light-element depletion due 
to transport by meridional circulation indicates that solar light-element abundance observations 
restrict the penetration to O.62i?0 (Charbonneau 2007); however this analysis does not necessarily 
suggest that the flow does penetrate that deep. Also vexing is the fact that different inversions, 
involving different helioseismic techniques such as ring-diagram or time-distance analysis recovers 
different profiles and widely varying peak meridional flow speeds (Giles et al. 1997; Braun & Fan 
1998; Gonzalez-Hernandez et al. 2006; Gizon & Rempel 2008). In our analysis, we chose to use 
the Gonzalez-Hernandez et al. data because at present, this provides the (relatively) deepest full 
inversion of the flow within the SCZ. Chou and Ladenkov (2005) reported time-distance diagrams 
reaching a depth of O.79i20 but have not yet reported a full inversion that could be used on our 
simulations. 

We point out that there is an important consequence of the presence of the flow speed maximum 
inside the convection zone - which is related to mass conservation: If the maximum poleward flow 
speed is found to be deeper inside the convection zone this would result in a stronger mass flux 
poleward, which needs to be balanced by a deeper counterflow subject to mass conservation; the 
density of the plasma increases rapidly as one goes deeper; e.g., the density at O.97i?0 is ten 
thousand times larger than at the surface. Although that is not achieved currently, our extensive 
efforts to fit the data point out that stronger constraints on the return flow may be achieved even 
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with data that does not necessarily go down to where this return flow is located, a fact that may be 
usefully utilized when better depth-dependent helioseismic data on meridional circulation becomes 
available. 

Although the depth of penetration of the circulation is an important constraint on the flow 
itself, our results indicate that the period of of the dynamo cycle does not in fact depend on this 
depth. Rather, our simulations point out that the period of the dynamo cycle is more sensitive to 
changes on the speed of the counterflow than changes anywhere else in the transport circuit, as this 
is where the dynamo loop originates. An accurate determination of the average meridional flow 
speed over this loop closing at the SCZ base is very important in the context of the field transport 
timescales. As shown by the analysis of Yeates, Nandy & Mackay (2008), the relative timescales 
of circulation and turbulent diffusion determines whether the dynamo operates in the advection or 
diffusion dominated regime - two regimes which have profoundly different flux transport dynamics 
and cycle memory (the latter may lead to predictability of future cycle amplitudes). Getting a firm 
handle on the average meridional flow speed is therefore very important and that is not currently 
achieved from the diverging helioseismic inversion results on the meridional flow. 

This suggests that a concerted effort using different helioseismic techniques on data for the 
meridional flow over at least a complete solar cycle (over the same period of time) may be neces- 
sary to generate a more coherent picture of the observational constraint on this flow profile. It's 
important to note that even though we used time averaged data, nothing prevents one from using 
the same methods to assimilate time dependent helioseismic data at different phases of the solar 
circle, allowing us to study the impact of time varying velocity flows on solar cycle properties and 
their predictability. 

On the encouraging side, our dynamo simulations show that it is relatively straightforward to 
use the available helioseismic data on the differential rotation (on which there is more consensus 
and agreement across various groups) within dynamo models. Also encouraging is the fact that 
the type of solar dynamo model presented here is able to handle the real helioseismic differential 
rotation profile and generate solar-like solutions. Moreover, as evident from our simulations, this 
dynamo model also generates plausible solar-like solutions over a wide range of meridional flow 
profiles, both deep and shallow, and with fast and slow peak flow speeds. This certainly bodes well 
for assimilating helioseismic data to construct better constrained solar dynamo models - building 
upon the techniques outlined here. 
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A. Numerical methods 

In order to use exponential propagation, we transform our system of Partial Differential Equa- 
tions(PDEs) in to a system of coupled Ordinary Differential Equations(ODEs) by discretizing the 
spatial operators using the following finite difference operators: 

For advective terms = + x(^))? where v is the velocity, we use a third order upwind 

scheme: 

dx \ g^(2A,+i + 3A-6A,_i + A,_2) iiv>0 ^ ^ ' ^ ' 

For diffusive terms — + x(^)^ 5 where r] is the diffusion coefficient, we use a second 
order space centered scheme: 



"d^ = (A^^^^-^ " + + ^^^^'^ ^^^^ 
For other first derivative terms = ff + xl^^)) we use a second order space centered scheme: 

dB 1 



dx 2A 



X 



(Ai-i - Ai+i) + 0{Ax') (A3) 



Here, x{^) corresponds to all the additional terms a PDE might have on the righthand side 
besides the term under discussion and Ai = A(xo + iAx)^i = 1,2, A^^ is our variable evaluated 
in a uniform grid of Nx elements separated by a distance Ax. 
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A.l. Exponential Propagation 

After discretization and inclusion of the boundary conditions we are left with an initial value 
problem of ordinary differential equations: 

^ = Fm)) (A4) 
U(to) = Uo (A5) 

where U is the solution vector in M^. Provided that the Jacobian dF{lJ{t)) exists and is continuous 
in the interval [tg, to + At], we can linearize F(U(to + At)) around the initial state to obtain 

= F(Uo) + aF(Uo)(U(to + At) - Uo) + i?(U(to + At)) (A6) 

where i?(U(to + At) are the residual high order terms. The solution to this equation can be written 
as 

AAt _ T 

U(to + At) = Uo + J + 0(At2) (A7) 

where A = dF{\Jo). Neglecting higher order terms leaves us with a scheme that is second order 
accurate in time and is an exact solution of the linear case. However there is a way of increasing the 
time accuracy of this method by following a generalization of Runge-Kutta methods for non-linear 
time- advancement operators proposed by Rosenbrock (1963). The combination of exponential 
propagation with Runge-Kutta methods was first proposed by Hochbruck and Lubich (1997) and 
then generalized by Hochbruck, Lubich and Selhofer (1998). In this work we use a fourth order 
algorithm which goes to the following intermediary steps to advance the solution vector between 
timesteps (U^) ^ U^+^)): 

ki = $QAtA^F(U^) 

k2 = $(AtA)F(U^) 
3 

'^3 0(^1 + ^2) 

o 

^3 = U'' + Atws 

ds = F{us) - FiJJ'') - AtAws 

ks = ^(^^AtA^ds 
Un+i = U- + At(^k2 + ^ks^ 

AAt _ T 

A = F(U^) , $ (AtA) = J 
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A.2. Krylov Approximation to the Exponential Operator 

Without any further approximation, this method is very expensive computationahy due to the 
need of continuously evaluating the matrix exponential. However, it is possible to make a good 
approximation by projecting the operator into a finite dimensional Krylov subspace 

Skt = span{Uo, AUq, A^Uq, A^-^Uq}. (A8) 

In order to do this we first compute an orthonormal basis for this subspace using the Arnoldi 
algorithm (Arnoldi 1951): 

1. i;i=Uo/|UoP 

2. For j = 1, m do 

• for z = 1, z compute hi^j — vj ^ Avj 

j 

• calculate w = Avj — hijVi 

• evaluate + 1, j) = \w\'^ 

• if hj^ij < e stop, else compute the next basis vector vj^i = w/hj^ij 

where e is the parameter that sets the error tolerance in this approximation. Once we have finished 
computing the algorithm, we have the relationship 

AV^ ^ VmH, (A9) 

where Vm = ['^i, •••,Vm] and H is a matrix whose elements are hij — vj ^ Avj. The validity of 
this approximation depends on the dimension of the Krylov subspace but numerical experiments 
have found that 15-30 Krylov vectors are usually enough (see Hochbruck and Lubich 1997; Tokman 
2001). Since {vi, ...,Vm} is an orthonormal basis of the Krylov subspace Sxr then V^Vm is a m x m 
identity matrix and VmV^ is a projector from onto Sxr- Using this projector we can find an 
approximation to any matrix vector multiplication by projecting them onto the Krylov subspace 
by calculating Ab ^ VmV^AVmV^b. After using Equation 



A9 



this becomes 



Ab ^ VmHV^b. (AlO) 

In fact we can approximate the action of any operator $(A), that can be expanded on a Taylor 
series, on the vector Uq using the Krylov subspace projection: 

$(A)Uo « Vm^H)V^Uo = |Uo| V^$(i?)ei (All) 



where ei — [1 ... 0] and we used vi = Uo/|Uo| and V^vi — ei. As we can see in Equation All 
the use of the Krylov approximation effectively reduces the size of the matrix operator; this makes 
the use of the exponential operator relatively inexpensive computationally. 



These two algorithms, in combination with a robust error control and an adaptative time-step 
mechanism strategy, form the core of the SD-Exp4 integrator (for more details see Hochbruck and 
Lubich 1997; Hochbruck, Lubich and Selhofer 1998 and Tokman 2001). 

In order to verify the performance standards of the SD-Exp4 code we ran case C of the dynamo 
benchmark by Jouve et al. (2008). This case is similar in nature to the simulations done in this 
work with the following differences: 

• The meridional flow profile has different radial and latitudinal dependence. 

• The poloidal source term has no quenching term and a different radial and latitudinal depen- 
dence. 

• The turbulent diffusivity profile consists of only one step and reaches a peak value of 10^^ 
crm? / s. 

• The analytic differential rotation has no cos^{9) dependence and uses a thinner tachocline. 



In order to compare the performance of the different codes in the benchmark study, two 
quantities are used: Cf^'^^ — a^'^^RQ/ricz which quantifies the minimum strength of the source- 
relative-to-dissipation that yields stable oscillations and uo — 2ti / [Trjcz] which quantifies the 
frequency of the magnetic cycle relative to the diffusion timescale. The dependence of this quantities 
is then plotted versus resolution for all the different codes (see Figure 11 of Jouve et al. 2008). In 
order to compare our code, we perform simulations with the same parameters and model ingredients 
as in case C of Jouve et al. (2008) and plot this two quantities versus resolution. We also plot 
the location of the mean found by the other codes and a shaded area encompassing one standard 



deviation around the mean. As can be seen in Figure 14, we find lower values of Cf^^ than the 
values obtained in Jouve et al. 2008 (1.86 as opposed to an average 2.46) for oscillatory solutions; 
this may be due to a lower amount of numerical diffusion in our code. Moreover, we find values of 
u in very good agreement with those found in Jouve et al. 2008 (540 as opposed to an average of 
538.2). We also find that for our code this quantities are less sensitive to resolution when compared 
to other codes used in the benchmark. This is likely due to the fact that we try to avoid numerical 
differentiation whenever we can do it analytically. 

As a final remark. Table 4 of Jouve et al. 2008 specifies the different time steps that are used 
by the different codes, which range from values of 10~^ to 10~^ code time units. Thanks to the use 
of Krylov approximations we are able to use time steps as large as 10~^ in code time units while 
solving Case C However, this does not translate to a performance improvement of several orders 
of magnitude (due to the added computations in the Krylov approximation). In order to quantify 
the relative performance improvement, we also made comparisons with the Surya code, which has 
been studied extensively in different contexts (e.g. Nandy & Choudhuri 2002, Chatterjee, Nandy 
& Choudhuri 2004, Choudhuri, Chatterjee & Jiang 2007), and is made available to the public on 
request. In comparison to the Surya code, we find that the SD-Exp4 code achieves a performance 
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improvement that reduces runtime from a half to a tenth of the total runtime depending on the 
particularities of the simulation. 
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Fig. 1. — Latitudinal(a) and radial(b) dependence of the poloidal source as quantified by the 
dynamo a-term. (c) Magnetic quenching of the poloidal source term. 



IVIagnetic Diffusivity Profile 




Fig. 2. — Radial dependence of the turbulent magnetic diffusivity. 
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Splines Interpolation of RLS Inversion 



Analytical Profile 
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Fig. 3. — (a) Spline interpolation of the RLS inversion, (b) Analytical profile of Charbonneau et 
al. (1999). (c) Differential rotation composite used in our simulations, (d) weighting function used 
to create a composite between the RLS inversion and the analytical profile of Charbonneau et al. 
1999. For all figures the red denotes the highest and blue the lowest value and the units are nHz 
with the exception of the weighting function. 
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Data Minus Composite Data Minus Analytical Profile 




0.2 0.4 0.6 0.8 1 (nHz) 0.2 0.4 0.6 0.8 1 (nHz) 



r/R r/R 

s s 

(a) (5) 

Fig. 4. — (a) Residual of subtracting the composite used in this work to the RLS inversion, (b) 
Residual of subtracting the analytical profile commonly used by the community to the RLS in- 
version. Red color corresponds to the hightest value and blue to the lowest. Graphs in units of 
nHz 
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Fig. 5. — Latitudinal velocity as a function of 9 used by van Ballegooijen and Choudhuri 1988 for 
different depths. Notice that the curves differ from each other only on their amplitude. 
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Meridional Flow Speed at Different Dephts Density Weighted Meridional Flow Speed at Different Dephts 
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Fig. 6. — (a) Measured meridional flow as a function of latitude at different depths (Courtesy 
Dr. Irene Gonzalez-Hernandez), each combination of colors and markers corresponds to a different 
depth ranging from O.97i?0 to Rq. (b) Meridional flow after being weighted using solar density. 
If we sum all data points at each latitude we obtain the average velocity, (c) Normalized average 
velocity and analytical fit. 
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Meridional Flow Speed for Different Latitudes Vth/G(th) for Different Latitudes 




(c) 

Fig. 7. — (a) Measured meridional flow as a function of radius at different latitudes (Courtesy 
Dr. Irene Gonzalez-Hernandez), each combination of colors and markers corresponds to a different 
latitude varying from —52.5 to 52.5. (b) Meridional flow after removing the latitudinal dependence, 
i. e. vo/G{6). The horizontal line with zero latitudinal velocity corresponds to the Equator, (c) 
Radial dependence of the latitudinally averaged meridional flow for our helioseismic data is depicted 
as large black dots. Other curves correspond to the radial dependence of the meridional flow profiles 
used in our simulations and solar density: Set 1 (black dotted) Rp — O.64i?0, Vq — 12m/ s] Set 2 
(magenta solid) Rp — O.64i?0, Vq — 22m/s] Set 3 (green dash-dot) Rp — O.71i?0, Vq — 12m/s and 
Set 4 (blue dashed line) Rp — O.71i?0, Vq — 22m/ s. The solar density taken from the solar Model S 
(Christensen-Dalsgaard et al. 1996) is depicted as a sohd red hne. The left-vertical axis is in units 
of velocity and the right-vertical is in units of density. 
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Fig. 8. — Butterfly diagram of the toroidal fleld at the bottom of the convection zone (color) with 
radial fleld at the surface (contours) superimposed on it. Each row corresponds to one of the 
different meridional circulation sets. The left column corresponds to runs using the helioseismic 
composite and the right one to runs using the analytical proflle. 
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Radial Composite Minus Analytical Shear Latitudinal Composite Minus Analytical Shear 
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Fig. 9. — (a) Residual after subtracting the radial shear of the analytical profile commonly used 
by the community from the radial shear of our composite data, (b) Residual of subtracting the 
latitudinal shear of the analytical profile commonly used by the community from the latitudinal 
shear of our composite data. 
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Fig. 10. — Snapshots of the shear source terms and the magnetic field over half a dynamo cycle (a 
sunspot cycle). Each row is advanced by an eight of the dynamo cycle (a quarter of the sunspot 
cycle) i.e., from top to bottom t = 0,T/8,T/4 and 3t/8. The solution corresponds to the composite 
differential rotation and meridional flow Set 1 (deepest penetration with a peak flow of 12 m/s) 
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Fig. 11. — Snapshots of the shear source terms and the magnetic field over half a dynamo cycle (a 
sunspot cycle). Each row is advanced by an eight of the dynamo cycle (a quarter of the sunspot 
cycle) i.e., from top to bottom t = 0,T/8,T/4 and 3t/8. The solution corresponds to the composite 
differential rotation and meridional flow Set 4 (shallowest penetration with a peak flow of 22 m/s) 
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Fig. 12. — Butterfly diagram of the torodial field at the bottom of the convection zone (Color) 
with radial field at the surface (Contours) superimposed on it when using a low diffusivity in 
the convection zone {rfcz — lO^^cm^/s). Each row corresponds to one of the different meridional 
circulation sets. The left column corresponds to runs using the helioseismology composite and the 
right one to runs using the analytical profile. 
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Fig. 13. — Snapshots of the magnetic field over half a dynamo cycle (a sunspot cycle) when using a 
low super-granular diffusivity {rjcz = 10^^ cm? /s). Each row is advanced by an eight of the dynamo 
cycle (a quarter of the sunspot cycle) i.e., from top to bottom t = 0,t/8,t/4 and 3t/8. The 
solutions correspond to the meridional flow Set 2 (deepest penetration with a peak flow of 12 m/s) 
and analytic differential rotation (left) and composite data (right) 
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Fig. 14. — Simulation results (dots) of running case C of the dynamo benchmark by Jouve et al. 
(2008). The top figure shows how the minimum value of Ca = aoRQ/rjcz foi" which the dynamo has 
stable oscillations depends on resolution. The bottom figure shows how u — 2tiR^/ {Trjcz) depends 
on resolution. In order to make the plot comparable to the benchmark, we include the mean value 
found by the different codes and a shaded area corresponding to one standard deviation around 
this mean. This run was made using the same meridional flow and poloidal source profiles as in 
the benchmark study (for details see Jouve et al. 2008). 



